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ABSTRACT 

The presence of a dark matter core in the central kiloparsec of many dwarf galaxies 
has been a long standing problem in galaxy formation theories based on the standard 
cold dark matter paradigm. Recent simulations, based on Smooth Particle Hydrody- 
namics and rather strong feedback recipes have shown that it was indeed possible to 
form extended dark matter cores using baryonic processes related to a more realistic 
treatment of the interstellar medium. Using adaptive mesh refinement, together with 
a new, stronger supernovae feedback scheme that we have recently implemented in the 
RAMSES code, we show that it is also possible to form a prominent dark matter core 
within the well-controlled framework of an isolated, initially cuspy, 10 billion solar 
masses dark matter halo. Although our numerical experiment is idealized, it allows 
a clean and unambiguous identification of the dark matter core formation process. 
Our dark matter inner profile is well fitted by a pseudo-isothermal profile with a core 
radius of 800 pc. The core formation mechanism is consistent with the one proposed 
recently by Pontzen & Governato. We highlight two key observational predictions of 
all simulations that find cusp-core transformations: (i) a bursty star formation history 
(SFH) with peak to trough ratio of 5 to 10 and a duty cycle comparable to the local 
dynamical time; and (ii) a stellar distribution that is hot with v/cr ~ 1. We compare 
the observational properties of our model galaxy with recent measurements of the iso- 
lated dwarf WLM. We show that the spatial and kinematical distribution of stars and 
HI gas are in striking agreement with observations, supporting the fundamental role 
played by stellar feedback in shaping both the stellar and dark matter distribution. 

Key words: galaxies: formation, dwarf - cosmology: dark matter - methods: numer- 
ical 



1 INTRODUCTION 

Dwarf galaxies, although very numerous and common in our 
present day universe, are also very faint and difficult to ob- 
serve. Nevertheless, it is now established that star formation 
proceeds at a very inefficient rate in dwarf galaxies, making 
them ideal laboratories to study the spatial distribution of 
their parent dark matter halo. Indeed, if dwarf galaxy are 
dark matter dominated, a stellar kinematic analysis gives 
direct constraints on the dark matter mass distribution. Al- 
though theoretical predictions of pure N-body models favor 
the formation of a cusp in the inner region of dark matter ha- 
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los ( Moore|1994 Navarro et al.|1 997), the observed rotation 
curve of dwarf and low surface brightness galaxies was shown 
to be more consistent with a shallower profile, even a con- 



stant density core (|de Blok fc McGaugh 1997 ; |de Blok et al. 



200Tl|de Blok fe Bosma|2002| |Kuzio de Naray et al.||2008 ) 
This apparent disagreement led to the so-called "cusp-core" 
problem, a serious challenge for the currently favored Cold 



Dark Matt er (CDM) paradigm ( |Moore|[T994| |Flores fc P"ri 
mack|1 994). Many solutions to this problem have been pro- 
posed in the recent years - for example a warm dark matter 
particle ( Kuzio de Naray et al.|2010 Villaescusa- Navarro & 
Dalai 2011 Maccio et al.|2012| or self-interacting dark mat- 
ter flSpergel fe Steinhardt|2000| |Rocha et al.|2012) - but the 
only explanation consistent with CDM relies on the effect of 



2 R. Teyssier et al 



baryonic physics, and more precisely on stellar feedback, to 
modify substantially the dark matter distribution. 

Mashche nko et alT| ( |2008| were the first to find a cusp- 
core transformation effect in a cosmological simulation, 
modeling a dwarf galaxy of mass ~ 10 9 M down to red- 
shift z ~ 5. This has been seen again in more recent work at 



higher masses and reaching lower redshifts ( Governato et al. 



2010||Macci6 et al.| [2Q12; Ma rtizzi et al.||2012| ). A long his- 
tory of work has looked at the possibility of baryonic physics 
generating such a transformation: for instance Navarro et al. 
(1996) showed that impulsive mass loss leads to irreversible 
expansion of orbits near the centre, Read &; Gilmore (2005) 
suggested that repeated epochs of outflows could produce 
a strong enough effect to generate cores and Mashchenko 
et al. ( 2006| pointed out that internal bulk motions of gas 
(not necessarily outflows) could have the same effect. Fi- 
nally, |Pontzen fc Governato] ( |2012| ) (hereafter PG12) pro- 
duced an analytic model of impulsive heating which was able 
to make quantitative predictions for the rate of cusp flatten- 
ing, agreeing with the simulated results of Govern ato et al.| 
(2010). This demonstrated that impulsive gas motions were 
the dominant cause of the simulated cusp-core transforma- 
tions as opposed to, for example, heating due to dynamical 
friction (e.g. |E1-Zant et al.|2001| |Goerdt etaT|2010 ). 

Although cosmological simulations provide realistic en- 
vironments and mass accretion histories for the model galax- 
ies, they are often challenging to analyze because of their ge- 
ometrical and evolutionary complexity. Moreover, since they 
are time consuming, the force resolution is currently limited 
to ~ 80 pc, which is only one tenth of the measured core ra- 



dius ( Governato et aL| 2010). It is therefore quite important 
to perform additional simulations with both a simpler set 
up and a higher resolution, which is precisely the goal of our 
paper. In Mashchenko et al. (2006), a similar approach was 



proposed with an isolated NFW halo stirred by 3 arbitrarily 
moving gas clumps. Although the proposed set-up was quite 
simplistic, it provided a well controlled experiment of cusp- 
to-core evolution with a better force resolution of 40 pc. Very 



recently, Cloet-Osselaer et al.| ( |2012| ) have performed a series 
of isolated NFW halo simulations with a more realistic treat- 
ment of star formation, stellar feedback and the induced gas 
motions based on the SPH technique. This work, following 
up on a series of papers reproducing many properties of the 



gas and stellar distribution in dwarf galaxies ( jValcke et al 
2008; Sch royen et al.|2011 ), also confirmed the formation of 
a dark matter core. 

In this paper, we will also perform idealized simulations 
of an isolated NFW halo, with a more realistic treatment of 
the ISM physics, and, for the first time, with the Adaptive 



Mesh Refinement (AMR) code RAMSES (Teyssier 2002}. It 
is indeed important to verify that core formation can also be 
recovered using a different type of code than the one used so 
far in all the previously cited papers, namely SPH in two dif- 



ferent implementations, GADGET by Springel (2005) and 



GASOLINE by |Wadsley et al.| ( |2004|) As shown by |Agertz 



et al. (2007) and Wadsley et al. (2008), both code types suf- 



fer from different systematic effects that might affect the 
numerical solution, especially when trying to resolve the 



clumpy ISM (see also |Read et al. (2010) and Read & Hay- 
ne!d|([20T2|)). 



The paper is organized as follows: in the first section, 
we describe the set up of our numerical experiment, with 



3 different simulations that we would like to compare. We 
then present and discuss our new stellar feedback imple- 
mentation in the RAMSES code, which allows us to have a 
much stronger dynamical effect on the surrounding gas than 
we were previously able to achieve, without invoking extra 
sources of energy. In the third section, we present our results, 
with emphasis put on the dark matter distribution and how 
the formation of the core correlates with strong potential 
fluctuations triggered by powerful gas expulsion phases. Fi- 
nally, we discuss and interpret our results in the light of the 
recently proposed mechanism by PG12, and confront our 
findings with observations. 



2 SIMULATION SET-UP 

We consider an isolated halo in hydrostatic equilibrium 



with / g , 



15 %. Both gas and dark matter follows the 



same NFW density profile, with concentration parameter 
c = 10. The halo circular velocity is chosen to be V200 = 
35 km/s, corresponding to the virial radius R200 = 50 kpc 
(or 35 kpc/h) and the virial mass M 20 o = 1.4 x 10 10 M© (or 
10 10 Mq/K). The halo is truncated at 112.5 kpc, so that the 
total enclosed mass is Mtot = 2 x 10 10 M Q . Following now 
rather standard prescriptions, such as in |Kaufmann et al7| 
( |2007| and pubois fe Teyssier | ( [20081 , we initialized the gas 
temperature by solving the hydrostatic equilibrium equation 
(see Eq. 1 in Kaufmann et al^]|2007[ for example), and we 
set the gaseous halo in slow rotation around the z-axis, us- 
ing the average angular momentum profile computed from 
cosmological simulation ( Bullock et al.|2001 ) and a spin pa- 
rameter A = 0.04. For more details on the set-up, the reader 
can refer to |Dubois fe Teyssier| (2008 ) . We use 1 million dark 
matter particles to sample the phase space of our dark mat- 
ter halo. Initial positions and velocities were computed us- 
ing the density-potential pair approach of Kazantzidis et al.| 
fl2004|), later refined by|Read et al.| (|2006p. To test the sta- 



bility of our gas-dark matter equilibrium system, we ran a 
first simulation with only adiabatic gas dynamics, without 
cooling nor star formation. As can be seen in Figure [I] our 
dark matter halo is stable during 2 Gyr of evolution and does 
not deviate by more than 10% from the initial NFW distri- 
bution, except within 60 pc from the center (corresponding 
roughly to twice the cell size), where the deviation can be 
as high as a factor of 2. 

All the simulations presented here have been run using 
the RAMSES cod e ( |Teyssier|2002| |Teyssier et al.|2006| [Fro 



|mang et al.||2006| . We used a quasi-Lagrangian refinement 
strategy to build the initial AMR grid, as well as refining and 
de-refining cells during the course of the simulation: each cell 
is individually refined if it contains more than 8 dark matter 
particles or if it contains a baryonic mass (gas mass + star 
particle mass) larger than 8 x m res , where m res = 1500 M©. 
Note that our gas mass resolution corresponds initially to 
2 million gas resolution elements across the whole gaseous 
halo. The box size being set to L = 300 kpc (in our ter- 
minology this corresponds to the first AMR level i — 1), 
we used isolated boundary conditions for the Poisson solver 
and zero-gradient boundary conditions for the hydro solver. 
For the latter, we used the HLLC Riemann solver and the 



MinMod slope limiter (see Fromang et al. ( 2006 ) for related 
details) . 
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Without cooling and star formation, we do not need to 
specify a maximum level of resolution: we just let the code 
refine the grid until it runs out of mass in the smallest AMR 
cell. In our adiabatic simulation, we have indeed reached 
"only" level t — 13 at the very centre of the halo, corre- 
sponding to a cell size Ax = 36 pc. If cooling and star forma- 
tion are activated, gas dynamics becomes strongly dissipa- 
tive and nothing can prevent the collapse. It is therefore cru- 
cial to set a maximum level of refinement, or a minimum grid 
resolution. In this paper, we use Umax 

— 14 or Ax — 18 pc, 
which is only a factor of 2 better than in the dark mat- 
ter only case, therefore avoiding two-body relaxation in the 
dark matter component. Our limited spatial resolution re- 
quires a careful treatment of our thermal model, in order to 
avoid numerical difficulties. As is now customary in galaxy 
and star formation simulations, we use an artificial pressure 
floor designed to enforce that the effective Jeans length is 
equal to 4 computational cells ( |Truelove et al.p997 ). This 
extra-pressure, equal to Pj ~ lQGAx p , is added to the 
thermal pressure in the Euler equation, gas cooling being 
applied only to the thermal internal energy. We use stan- 
dard H and He cooling ( Kat z et al.|1996|, with an additional 
contri bution from metals based on the [Sutherland & Dopita 
( |l993| ) model above 10 4 K and metal fine-structure cooling 
below 10 4 K as in Rosen & Bregman (1995). Metallicity is 



modeled as a unique, passively advected quantity, noted Z, 
representing the mass fraction of metals, and seeded by in- 
dividual supernovae event with a yield y = 0.1. The initial 
metallicity in the halo is set to Z — 10 -3 Z©, enabling from 
the start cooling down to very low temperatures. With these 
ingredients and for our spatial resolution of 18 pc, gas will 
cool efficiently down to Tj c± 600 K and reach a density 
of rtj ~ 60 H/cc, before hitting the Jeans- length related 
pressure floor (the index J standing for Jeans-length related 
quantities) . This is a conservative estimate of the maximum 
gas density we can reach before being affected by finite res- 
olution effects. 

Star formation (SF) is modeled using a Schmidt law 
( Schmidt||1959| ), for which the SF rate is given by 



/9* — 6* 



if />gas > P* 



(i) 



The SF efficiency per free-fall time t& is usually chosen close 



to a few percent (Krumholz & Tan 2007); we used here 
e* = 0.01. The density threshold for SF, above which gas is 
considered to be dense enough to be eligible to form stars, 
has to be chosen carefully. In order to mimic as closely as 
possible realistic SF in dense molecular clouds, we want to 
choose this threshold as high as possible. On the other hand, 
as discussed previously, the gas density distribution in the 
simulation will be affected by finite resolution effect close to 
the Jeans density, in our case nj c± 60 H/cc. It is therefore 
quite natural to choose the SF threshold density close to the 
Jeans density. In the present paper, we conservatively picked 
n* = nj/3 = 20 H/cc. These recipes have been all imple- 
mented in the RAMSES code and used for many galaxy 



formation studies in the recent years ( 


Agertz et al. 2009; 


Teyssier et al.||2010| |Bournaud et al.| 


2011| |Agertz et al.| 


2011 


Scannapieco et al.|2012). We will now present in more 



details the main ingredient of the present study, namely our 
new implementation of stellar feedback. 



100.00 



10.00 



1.00 



0.10 



0.01 




t=0.0 Gyr with NFW 
t=0.5 Gyr 
t=1.0 Gyr 
t=1.6 Gyr 



0.01 



0.10 



1.00 



10.00 



(kpc) 



Figure 1. Evolution of the dark matter density profile over the 
2 Gyr of evolution for the control run with adiabatic (no cooling 
and no star formation) hydrodynamics. We see that the dark 
matter halo density profile remains stable (with at most 10% 
deviations), except within the central 60 pc, which corresponds 
to twice the cell size in this case, in which we see less than a factor 
of two deviations from the initial profile. 



3 STELLAR FEEDBACK IMPLEMENTATION 

The last but probably most important ingredient in galaxy 
formation simulations is stellar feedback. Stars, especially 
young stars, interact quite strongly with their environment. 
Among many possible physical processes, HII regions, stellar 
winds, type II and type la supernovae explosions have been 
the most popular. These various mechanisms are believed 



to be responsible for SF regulation in dwarf galaxies (Dekel 
&; Silk|p"986 ), and they contribute collectively to roughly a 



couple of 10 51 erg of energy injection into the surrounding 
ISM, for an average 10 Mq massive star and for a standard 



stellar initial mass function (IMF) ( |Piontek Sz Steinmetz 
|2011| ). Since all these processes occur at very small, sub- 
pc scales, we have to rely on rather approximate numerical 
implementations. The simplest approach is clearly to inject 
directly thermal energy into each gas cell containing young 
stars. It turns out that feedback is very inefficient in that 
case because it occurs mostly in dense gas where cooling 
immediately radiates away the added thermal energy. It is 
usually considered that this strong cooling is a spurious nu- 
merical effect ( |Ceverino &; Klypm" 2009). But supernovae 
explosions, and more generally feedback from young stars, 
takes place in very dense environment (such as molecular 
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Figure 2. Maps of gas density (left, in logarithmic scale between 0.1 et 1000 H/cc), gas temperature (middle, in logarithmic scale 
between 300 and 10 6 K) and stellar column density (right, logarithmic scale in arbitrary units), seen face-on (lower plot) and side-on 
(upper plot), with (upper half) and without (lower half) feedback after 2 Gyr of evolution. The images are all 10 kpc across. 
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clouds) so that strong radiative effects are likely to take 
place anyway. More recently, radiative pressure from young 
stars has been also proposed as a new possible mechanism 
(Mu rray et al.pblO ) to drive momentum into the interstel- 
lar medium (ISM). This approach is very promising since 
gas momentum does not suffer directly from gas cooling. It 
is however quite difficult to model since it requires detailed 
radiative transfer calculation of the infrared radiation repro- 
cessed by dust. It was also shown recently that, for dwarf 
galaxies such as the ones studied here, super novae feedback 



likely dominates over other sources of feedback energy (Hop 
kins et al.|2012| ) : radiative effects become inefficient due to 
the lower column densities. 

Physical processes associated to stellar feedback are 
probably even more complex than the ones quoted here so 
far. For example, we know from X-ray observations of su- 
pernovae remnants than they are strongly magnetized and 
turbulent. Moreover, cosmic rays acceleration occurs very 
rapidly and the energy density stored in the relativistic com- 
ponent is large enough to affect significantly the dynamics 
of the propagating shock wave ( Ellis on et al.||2004 Ferrand 
et al.|20To| . These non-thermal processes have much longer 
dissipation time-scales than the thermal component, so they 
can store feedback energy longer than the rapidly cooling 
gas and release it to the gaseous component more gradually. 
Modeling unresolved turbulence, magnetic field amplifica- 
tion and cosmic ray propagation and radiative losses in the 
framework of galaxy formation is at its infancy, although 
some recent progresses hav e been reported (|Hanasz et aT 
5006] |Jubelgas et al.| [20081 |Scannapieco fe Briiggen||2010| 
Uhlig et al.|2012) . We propose here a much simpler formal- 
ism that captures very crudely these various non-thermal 
processes and couple more efficiently the energy associated 
with stellar feedback to the gas component. 

The idea is to introduce a new variable for the energy 
density in these non-thermal components. We note this new 
variable e tU rb for simplicity, although it doesn't have to be 
associated to turbulence only, but also to cosmic rays and 
magnetic fields. The specific turbulent energy eturb is defined 
by e t urb = peturb- We model the time evolution of this non- 
thermal energy using the following simple equation 



De t 



Dt 



P^turb 
^diss 



(2) 



where we have only 2 terms, the no n- thermal energy source 
Ei n j due to stellar feedback and the energy dissipation mod- 
eled classically as a damping term with dissipation time- 
scale tdiss- The time scale depends on the exact underlying 
dissipation mechanism, such as radiative losses for cosmic 
rays ( |EnBlin et al.|[2007| or eddy turn-over time for turbu- 
lence. In what follows, we just need to accumulate enough 
feedback energy to drive strong outflows and to regulate our 
SF efficiency across the galactic disk. In this paper, we used 
a fixed dissipation time scale 



tdi. 



10 Myr 



(3) 



comparable to the typical molecular cloud life time, as ex- 
pected from their observed internal SF efficiency (IWilliamsl 



& McKee 1997). Other options are possible, such as, for 
example, tdiss — Aj/cr tU rb, where Aj could be the typical 
Jeans length or, equivalent ly, the typical cell size, or even 
more exotic, the cosmic rays mean free path or magnetic 



dissipation length scale. The turbulence velocity dispersion 
crturb is defined here as 



_ 1 2 

tturb — ^ ^turb- 



We parametrized the non- thermal energy injection as 
E inj = /j*ry S Nl0 50 erg/M , 



(4) 



(5) 



which corresponds to a maximally efficient model for which 
10 51 erg per 10 M© massive star is injected into non-thermal 
energy. For a typical IMF, the mass fraction of stars going 
supernovae is roughly ?7sn = 0.1. Although these numbers 
are determined by stellar evolution models, we believe that 
they can be considered as free parameters. Indeed, both the 
IMF and stellar evolution models are poorly constrained at 
low metallicity and/or at high redshift, and more generally 
in extreme galactic environments ( Marks et al.|[2012 ). Note 
that this energy is injected in the gas cell only 10 Myr after 
the star particle has been created. 

The final step of our feedback model is to compute the 



non-thermal pressure using P t urb 
thermal pressure to get 

-ftot — -^thermal ~\~ -Pturb- 



• pcr turh and add it to the 



(6) 



This obviously requires to modify the Euler equation by 
adding the non-thermal pressure term, and this could im- 
pact more or less deeply the hydrodynamics solver. One even 
simpler alternative is to actually add directly to the thermal 
energy the injected non-thermal energy as 

-D6thermal -A- 



Dt 



Pthe 



iV • v - n H A 



(7) 



and shut down gas cooling to mimic the contribution of 
the non-thermal pressure to the total pressure, eturb be- 
comes a completely independent variable (like a passive La- 
grangian tracer) used only to re-activate cooling when the 
non-thermal contribution becomes comparable or smaller 
than the thermal energy. In practice, we shut down cool- 
ing everywhere the turbulence velocity dispersion is large 
enough, using 



A = if cr t urb > 10 km/s. 



(8) 



Our final model does not require any modification to the 
hydro solver. It just requires to shut down cooling in each 
cell, if the velocity dispersion is above a chosen threshold 
(here 10 km/s). The other advantage is that it is qualita- 
tively very similar to the "delayed cooling" approach used 
in other recent studies (Stinson et al. 2006; Governato et al. 



2010 Agertz et al.||2011 ), although our current approach is 
motivated by non-thermal astrophysical processes for which 
a detailed modeling would be far beyond the scope of the 
present study. 



4 EVOLUTION OF THE DARK MATTER 
DENSITY PROFILE 

Now that we have presented our numerical set-up and pa- 
rameters, we now describe in more details our results. We 
have performed 3 different simulations, in order to highlight 
the effect of stellar feedback on the mass distribution within 
the halo. The first simulation, already presented in the previ- 
ous sections, was a pure adiabatic gas evolution of our initial 
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Figure 3. Star formation history in the runs without (left plot) and with (right plot) feedback. 



hydrostatic equilibrium halo. We checked that our initial set- 
up was indeed stable over 2 Gyr time, as it should be since 
no gas cooling was considered in this case (see Fig. [I). We 
would like to stress that this is a very important step in our 
methodology, since it demonstrates that any evolution in the 
dark matter density profile has to be related to the dissipa- 
tive nature of baryons, through gas cooling, star formation 
or feedback. The second simulation was run with (metal de- 
pendent) gas cooling and star formation. No stellar feedback 
was included. In this case, the gas loses pressure support and 
rains down towards the centre of the halo, quickly reaching 
the centrifugal barrier and sets up into a centrifugally sup- 
ported disc. As one can see in Figure [2] the final disk is very 
thin and fragments into dense gas clumps that form stars ac- 
tively. This gives rise to the formation of dense bound star 
clusters that survive for long times. We also see the forma- 
tion of a massive bulge in the centre of the galaxy, leading to 
an overall highly concentrated baryons distribution. The as- 
sociated SF history can be seen on Figure [3] it is on average 
very high, around 1 Mo/yr, with short bursts reaching 4 to 
6 M /yr, associated to the formation of dense gas clumps. 
Such a high SF rate is usually associated to massive galaxies 
at low redshift. This is quite unrealistic for dwarf galaxies we 
see today (Hopkins et al. 2002). The effect of this strongly 
dissipative evolution on the dark matter profile can be seen 
on Figure [4] After 1 Gyr, the dark matter distribution has 
been adiabatically contracted very significantly by baryons. 
The inner slope of the dark matter density profile is close 
to -2, and no core is visible. It is worth mentioning that al- 
though we have a very clumpy structure in the ISM and in 
the stellar distribution, it does not trigger the formation of 



a dark matter core in our case: the mechanism proposed by 
Mashchenko et al.| ( [2006] ) (see also |El-Zant et al.|200l) does 
not work here, probably because our clumps are not massive 
enough. 

We now move to our final run with gas cooling, SF and 
stellar feedback. The evolution of the star forming disk is 
dramatically different from the "no feedback" run. We see 
in Figure [2] that the final gas distribution shows a very thick, 
turbulent disk, with strong outflows made of shredded clouds 
and filaments. The face-on view reveals that many gas clouds 
form in the outskirts of the disk, while the central region has 
been evacuated by stellar feedback, giving rise to the wind. 
The temperature map illustrates nicely the hot gas in the 
wind, segregating from the cold gas in the ISM. Star for- 
mation still proceeds within dense clouds, but these are not 
long-lived anymore. This is why we don't see any massive 
star clusters in the stellar surface density map. Only a few 
managed to survive. This is one of the key qualitative fea- 
tures of our stellar feedback implementation: gravitational 
instability and shock compression trigger the formation of 
star forming clouds, which are then quickly disrupted by 
stellar feedback, recycling the unused gas into the ISM and 
giving rise to the galactic wind. 

The SF rate plotted in Figure [3] is one order of mag- 
nitude lower in average than the "no feedback" case. It ex- 
hibits strong bursts followed by quiescent phases. When gas 
cools down and sinks towards the central region, SF rises 
sharply and triggers a starburst. Stellar feedback then re- 
moves the gas into the hot wind, leaving the central kpc 
almost devoid of gas. This explains the very low star forma- 
tion episodes. Gas then rains back down from the corona and 



Cusp- core transformations in dwarf galaxies 7 



100.00 



10.00 



55 



1.00 



0.10 



0.01 




0.01 



0.10 



1.00 



10.00 



(kpc) 



100.00 



10.00 



1.00 



0.10 



0.01 




0.01 



t=0.0 Gyr with NFW 
t=0.5 Gyr r core =500 pc 
t=1.0 Gyr r core =700 pc 
t=1.6 Gyr r core =800 pc 



0.10 



1.00 



10.00 



(kpc) 



Figure 4. Evolution of the dark matter density profile over the 
2 Gyr of evolution for the control run with cooling, star forma- 
tion but no feedback. The dark matter halo has been strongly 
adiabatically contracted. 



Figure 5. Evolution of the dark matter density profile over the 2 
Gyr of evolution for the control run with cooling, star formation 
and stellar feedback. We see the formation of a large core. We 
also show for comparison the analytical fit (dashed line) based on 
a pseudo-isothermal profile (see text for details) 



triggers a new SF episode. These cycles are clearly visible 
in the SF history. They are at the origin of strong poten- 
tial fluctuations due to massive periodic gas outflows and 
inflows. 

In Figure [5] we plot the evolution of the dark mat- 
ter density profile. As a technical side note, we would like 
to stress that these profiles have been computed using the 
highest dark matter density peak as centre, defined using 
the "shrinking sphere" algorithm. This minimizes spurious 
features in the density profiles due to poor centering. We 
clearly see in the "feedback" run that a large dark matter 
core develops. We tried to fit the resulting profile using a 
"pseudo-isothermal" profile defined as 



p oc 



1 + (r/r CO re) 2 



(9) 



This analytical shape has traditionally been used in rota- 
tion curve fitting and has proven to work extremely well in 
nearby dwarf galaxies ( Begeman et al.|1991| Kuzio de Naray 
et al. 2006, 2008). As can be seen in Figure [5j this is also 
true in our simulation, for which the pseudo-isothermal fit is 
remarkably good within the central 3-4 kpc. We have mea- 
sured core radii ranging from 500 pc at early time (0.5 Gyr) 
to 800 pc at late time (1.6 Gyr). We do see a systematic 
increase of the core radius with time. We have also tried to 
fit our dark matter profile with other popular cored den- 



sity distributions (Burkert 1995, for example) but with less 
success. 

To quantify even more the time evolution of the dark 
matter distribution, we have fitted the density profile be- 
tween 200 pc and 800 pc with a single power law. This is 
obviously a poor fit to our cored distribution, but this cap- 
tures the essence of the dark matter flattening in the centre. 
At t — 0, we measure for the slope a — —1.2, which is the 
average value of the NFW profile between the 2 chosen radii. 
At the beginning, gas cools down and adiabatic contraction 
of the halo can be measured as a decrease of the inner slope 
to a = —1.8. But stellar feedback processes very quickly de- 
velop and trigger strong potential fluctuations, leading to a 
gradual flattening of the slope. We end up after 2 Gyr of 
evolution with a slope between a = —0.3 and a = —0.5, 
in excellent agreement with the measurement reported by 
PG12 in the context of a lower resolution cosmological sim- 
ulation. 

As proposed by PG12, this strongly suggests that the 
core formation mechanism is related to violent gas outflows 
triggered by a succession of starbursts. In order to validate 
this idea in our simulation, we have computed the time evo- 
lution of the enclosed gas mass within spheres of increas- 
ingly large radii, from 200 pc for the smallest to 1600 pc for 
the largest. The results are shown in Figure 7 and illustrate 
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that gas is repeatedly removed from the central region of the 
galaxy, giving rise to strong coherent potential fluctuations. 
As more gas cools down from the halo, the SF bursts become 
more violent and large mass fluctuations propagate to larger 
radii. This explains why we see a systematic increase in the 
dark matter core size as a function of time. At the end of the 
simulation, only scales up to ^800 pc are affected by these 
strong fluctuations, corroborating nicely the measured core 
size of 800 pc. 

While these results are qualitatively consistent with the 
analysis of PG12, there is a difference in the quantitative de- 
tails. In our case, many potential fluctuations are not only 
established in under a dynamical time (as in PG12) but also 
are erased before a full dynamical time has elapsed (unlike 
in PG12). The effect of these short-time transients is ig- 
nored by the PG12 toy model, since it starts by sampling 
the potential only once per dynamical time. Applying the 
PG12 algorithm directly therefore leads to incorrect predic- 
tions for our present simulation. We verified, however, that 
consistent predictions can be made by changing the sam- 
pling method to pick out, at each dynamical time step, the 
most extreme potential jump that occurred within the sim- 
ulation. While this is a heuristic method, we believe it shows 
that the essential ingredient of impulsive potential changes 
is at play in our new simulations. It is in principle possi- 
ble to extend the PG12 work to track phases of particles, 
thus giving a complete description of the impulsive changes 
which does not require the problematic quantization of time. 
However this approach would lose the attractive simplicity 
of the PG12 model and is certainly beyond the scope of the 
present investigation. 

We should also emphasis that the core forms after only 
0.5 Gyr in our reference feedback run, when only 10% of the 
stars have formed. We have indeed measured a total stellar 
mass of M* ~ 5 x 10 7 at that time. This means that 
the final total energy liberated in supernovae explosions is 
not the limiting factor. This is in fact consistent with the 
estimates of energy required to form cores in Peharrubia 
et al. (2012). The question of what limits star formation 



requires full cosmological runs to tackle further, since our 
"monolithic collapse" models do not reflect the interplay 
between gas inflow, outflow and feedback that presumably 
regulate star formation as of function of redshift in the real 
universe. 



5 OBSERVATIONAL GALAXY PROPERTIES 

We have shown that cusp-core transformations are expected 
if there are strong (order unity) potential fluctuations on 
dynamical timescales. Such fluctuations are non-adiabatic, 
non-reversible, and heat both the dark matter and stars. We 
discuss here the observational consequences of such fluctu- 
ations; these can then be used to constrain or rule out the 
presence of such fluctuations in real systems. 



5.1 The stellar distribution 



Read & Gilmore (2005) found from their toy models that, 



when cusp-core transformations occur, the resulting surface 
brightness of the stars is well-approximated by an exponen- 
tial over many scale lengths, sometimes with a break radius 
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Figure 6. Time evolution of the slope of the dark matter density 
profile measured between 200 and 800 pc for the simulation with 
feedback. 



at large radii and a core at small radii. This is exactly what 
we see in the present simulation: in Figure [8] we have plot- 
ted the stellar surface density with and without feedback. 
In the latter case, we see a prominent bulge with high Sersic 
index and a weak exponential disk extending up to 6 kpc, 
while in the feedback case, we see a smaller exponential disk, 
with no bulge in the center and a clear signature of a stellar 
core within the central 500 pc. We have fitted this stellar 
distribution with a pure exponential disk with scale length 
Rd — 1.1 kpc (see Figure [8]). In our simulation without feed- 
back, the stars form in a razor-thin, rotationally supported 
disc. 

Note that simulations of the formation and evolution 
of isolated dwarf galaxies generically end up with near- 
exponential surface brightness profiles, with rather thick gas 
and stellar disks (Springel &; Hernquist 2003; Valcke et al 



2008 



Dubois & Teyssier 2008). These simulations can be 



however divided into 2 main categories: quiescent feedback 
models, in which small scale effects are captured by an ef- 
fective equation of state, leading to a quiescent star forma- 
tion history where outflows are obtained through a quasi- 
stationary galactic wind (see e.g. |Springel &; Hernquist|2 003; 
Dubois & Teyssier 2008); and bursty models, in which the 



star formation history is not stationary and shows violent 
fluctuations (see e.g. |Stinson et al.|2007| [Valcke et al.|2008| 



Revaz et al.|2009 ) . The first category of simulations give rise 
to moderately thick disk, and very weak potential fluctua- 
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Figure 7. Time evolution of the total enclosed gas mass within 
spheres of radius 200 (blue), 400 (green), 800 (red) and 1600 
(black) pc for the simulation with feedback. 



Figure 8. Stellar surface density of the final dwarf galaxy with 
(red line) and without (black line) feedback. In blue is shown the 
exponential fit of the stellar disk with feedback with scale length 
Rd = 1.1 kpc. 



tions. The stellar disk is still well defined with an aspect 
ratio close to h/R ~ 0.1 (see Fig. 11 in Springel & Hern- 



quist||2003[ for example). Therefore, we should not expect 
any dark matter core in this case. The second type of simu- 
lations, however, gives much thicker disks, in some case al- 
most spheroidal galaxies. They show violently time- varying 
outflows, and they should lead to the formation of a dark 
matter core. Our present model, with strong stellar feed- 
back, belongs to the second category: the resulting stellar 
distribution is oblate and quite hot, with v/a ~ 1 and a 
scale height to scale length ratio close to 0.5 (see Figure [5] 
and Figure 10). 



Unfortunately, it can be difficult to conclusively test the 
above expectations. Exponential surface brightness profiles 
can form through a variety of physical processes not related 
to violent potential fluctuations (e.g. Lin & Pringle 1987). 
Similarly, hot stellar distributions with v/a ^ 1 need not 
be a smoking gun for cusp/core transformations. Galactic 
mergers (e.g. |Searle &; Zinn 1978; Re ad et al.|2006 ), and col 



lisionless heating from a strong tidal field (e.g. Mayer et al. 



2001 Lokas et al. 2011) both lead to stellar distributions 



with low v/a. 

We can avoid all of the above complications, how- 
ever, if we focus on isolated low mass dwarf galaxies. The 
best-studied system to date is WLM. Despite being ex- 
tremely tidally isolated, WLM does indeed have a hot, 
oblate spheroidal, stellar distribution (with v/a ~ 1 and 



h/Rd — 0.5) that is reasonably approximated by an expo- 
nential surface brightness fall-off ( Leaman et al.|2012| . Note 
that in this case too, a flattening on the stellar surface den- 
sity profile is seen in the centre, reminiscent of a stellar cor^ 
iLeaman et al. (2012) point out that it may be difficult to 



understand why the stellar distribution in WLM is so hot 
without recourse to strong and bursty stellar feedback. We 
can confirm that view here. In our model of an isolated dwarf 
galaxy, only the simulation with very strong feedback man- 
aged to form a hot thick-disc-like stellar distribution. This 
same simulation found a significant cusp-core transforma- 
tion in the underlying dark matter distribution. In order to 
push the comparison with WLM even further, we have an- 
alyzed the kinematic properties of our dwarf galaxy. These 
are shown on Figure [l0| One can see immediately that the 
stellar rotation curve is slowly rising, while the stellar veloc- 
ity dispersion is slowly declining. Both curves intersect at 
around 2 kpc, at a value close to 20 km/s. These features 
are in striking agreement with WLM data as exposed in |Lea-| 
|man et aT] ( |2012| . The gas kinematic properties are also very 
similar, as shown in Figure [l0| Using our kinematic analy- 
sis, we are in a good position to test one important method 
to derive the total mass profile in dwarf galaxies from their 



1 PG12 pointed out the heating processes applying to the DM 
should also apply to the stars. 
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Figure 9. Kinematic analysis of our dwarf galaxy without feed- 
back: circular velocities (black solid line), stellar tangential ve- 
locity (red solid line) and stellar tangential velocity dispersion 
(green solid line), compared to the gas tangential velocity (red 
dotted line) and gas tangential velocity dispersion (green dotted 
line). 



Figure 10. Kinematic analysis in the feedback case: circular ve- 
locities (black solid line), stellar tangential velocity (red solid line) 
and stellar tangential velocity dispersion (green solid line), com- 
pared to the gas tangential velocity (red dotted line) and gas 
tangential velocity dispersion (green dotted line). Also shown as 
the blue solid line is the predicted circular velocity curve based 
on the Asymmetric Drift (AD) approximation. 



kinematic properties, namely the Asymmetric Drift (AD) 
model. This method, based on the Jeans equation, follows 
a few reasonable simplifying assumptions to derive the rela- 
tion between the circular velocity and the velocity moments. 
Following |Hinz et al. ( |2001 ) and |Leaman et ah] (2012), we 



used v circ — v e + (7(9 (2r/Rd — 1) with Rd ~ 1.1 kpc, 
measured in our simulation. We see in Figure [To] that AD 
is overall a good approximation to recover the underlying 
mass profile, except perhaps in the very center where it is 
unde restimated. The tot al mass inferred from this analysis 
by |Leaman et ah] ( |2012| ) for WLM, M tot ~ 2 x 10 10 M 
is therefore accurate, and again very close to our simulated 
halo mass. 

Although our spatial and kinematic properties are in 
striking agreement with the relatively isolated dwarf WLM, 
the total stellar mass that we obtained in our simulation 
is too large by one order of magnitude. We have plotted 
in Figure [TT] the cumulative stellar mass profile in spheri- 
cal shells. One sees clearly that without feedback, almost all 
baryons are converted into stars after 1 Gyr, since we get 
M* c± 10 9 M©. With our strong stellar feedback model, we 
managed to reduce this number by one order of magnitude, 
down to M* c± 10 8 M©. This is quite an achievement, but it 
falls short by one order of magnitude to explain the stellar 



mass observed i n WLM, which has been measured by|Jack- 
son et al.M2007k to be M* ~ 1.1 x 10 7 M . We are therefore 



overproducing stars to a level comparable to most current 



galaxy formation simulations (Piontek & Steinmetz 2009 



Governato et al.]|2010| |Agertz et al.||20TT ) , when compared 
to individual galaxies or to an ensemble of galaxies using 
the abundance matching technique ( Guo et al.|2010 



et al. 2010), although recently iMunshi et al. (2012) argue 



Moster 



differently. Although solving this issue is beyond the scope 
of the present paper, we have a conceptually simple way to 
solve this problem, by lowering the star formation efficiency 
parameter e* by one order of magnitude, and in the same 
time, increasing the mass fraction of massive star going su- 
pernovae by also one order of magnitude. The first idea could 
be justified by the low metallicity we find in dwarf galaxies, 
leading to a inefficient regime of star formation, for which 
dust shielding is less efficient at promoting H 2 molecule for- 
mation ( |Krumhofz fc Dekel|201l] ). The second idea could be 
justified by recent observations of low metallicity star clus- 
ters in the Galaxy, which are consistent with a top-heavy 
IMF ( Marks et al.|2012 ). Using these two non-standard but 
plausible ingredients, we will straightforwardly obtain the 
same energy input from supernovae, and therefore the same 
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Figure 11. Cumulative stellar mass profile as a function of spher- 
ical radius for the fiducial star formation efficiency (1%) without 
(black line) and with (red line) stellar feedback at t ~ 1 Gyr. Also 
shown are the stellar mass profiles with a very low star formation 
efficiency (0.2%) and a top heavy initial mass function, without 
(green line) and with (blue line) stellar feedback. Only the latter 
model can match the observed stellar mass in the dwarf galaxy 
WLM. 



form most of its baryons into stars. With stellar feedback, 
however, we form only M* ~ 2 x 10 7 of stars after 1 Gyr, 
in much better agreement with the observed stellar mass in 
WLM. We have also checked (not shown here) that, in this 
last simulation, all of the galaxy properties discussed above 
are also recovered, namely a large dark matter core, thick 
stellar and gaseous disks and a bursty SF history, with large 
amount of gas moving in and out of the central kpc of the 
galaxy. This last model gives very encouraging results, but 
requires a significantly top-heavy IMF together with a very 
low star formation efficiency. We would like also to stress 
again (see Section 4) that star formation might be also reg- 
ulated or even stopped after only 0.5 Gyr of evolution, due to 
the effect of a more realistic cosmological accretion history, 
which is not captured here in this idealized set-up. 



5.2 Star formation histories 

A second key observational prediction of cusp-core transfor- 
mations is the star formation history (SFH). In the category 
of feedback models where cusp-core transformations are ex- 
pected (|Stinson et aL|2007| |Valcke et aL|2008| [Revaz et al 



2QQ9; |Cloet-Osselaer et al.|2012[ and this paper), the SFH is 
extremely bursty with peak to trough variations of 10 and 
a duty cycle of roughly one dynamical time. Although our 
average star formation rate of 0.02 to 0.2 M©/yr for a total 
gas mass of 10 9 Mq is consistent with observations of blue 



compact dwarf galaxies in the local universe (jHopkins et al. 



2002), the bursty nature of its time evolution needs to be 
tested observationally. There are two possibilities: (i) mea- 
sure the SFH for individual systems; and (ii) estimate the 
variance in the star formation statistically for a large pop- 
ulation of like-galaxies. The former is cleaner, but requires 
very high time resolution in the derived SFH; the latter has 
potentially high time resolution but relies on assumptions 
about the equivalence of large populations of galaxies. We 
discuss each in turn, next. 



hydrodynamical model, with however ten times less long- 
lived star formed. A different model based on early radiative 
stellar feedback has been recently proposed by |Brook et al.| 
( 2012| that could solve this very same issue, although the 
quantitative effect would not be as straightforward as the 
one proposed here. 

Since complex non-linear hydrodynamical processes are 
at play, it is not obvious that lowering e* (while increasing 
the feedback energy per stellar mass formed) will result in 
decreasing the total stellar mass formed. We therefore run 
2 additional simulations with e* = 0.2% and a top-heavy 
IMF so that 50% of the mass of a given stellar population 
(or star particle) will be in massive stars going supernovae. 
The first simulation was run without stellar feedback (only 
metal enrichment) and the second simulation was run with 
our new stellar feedback scheme and a top-heavy IMF. We 
have plotted in Figure [TT] the cumulative stellar mass pro- 
files of these 2 additional runs. Interestingly, although we 
have decreased by a factor of 5 the SFE, the run without 
feedback gives also rise to the formation of M* c± 10 9 Mq 
of stars, like in the fiducial case. This means that, without 
stellar feedback, the gaseous disk always manages to cool 
and contract enough to reach high gas densities and trans- 



5.2.1 Individual star formation histories 

The best SFHs come from the most nearby systems. These 
have resolved color magnitude diagrams (CMDs) that reach 
the oldest main sequence turn-offs, which is vital for cor- 



rectly recovering the intermediate age stars (e.g. iNoel et al. 



2009). More distant systems have SFHs derived via spectral 
energy distribution (SED) fitting of the unresolved stars that 
is degenerate and less accurate than CMD fitting (e.g. Zhang 
pFaLl[2012| ). 

The best-studied system to date is the nearby dwarf 
spheroidal galaxy Sculptor that orbits the Milky Way (e.g. 
Mateo|1998||Lux et al.|2010| ). Unfortunately, even for excel- 
lent data, CMD fitting has a temporal resolution poorer than 
-lGyr ( |Dolphin|[2002t , where the errors come from a mix 
of photometric uncertainty, spread in distance modulus, and 
errors in the stellar population evolution libraries \ Dolphin 



2012). de Boer et al. (2012) have recently attempted to im- 
prove this situation by simultaneously fitting the CMD and 
spectroscopic metallicities of red giant branch stars. They 
explicitly consider how well their methodology could recover 
a bursty star formation history. Using simulated data, they 
show that their method would recover a smooth continuous 
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star formation history (as observed for Sculptor) from an 
input bursty history. Thus, the very best data remain in- 
conclusive. It is interesting to note, however, that the very 
latest dynamical models for Sculptor appear to favor a cored 
rather than cusped mass distribution ( Battaglia et al.| 2QQ8 



Wal ker fc Penarru bia 2011 |Amorisco Evans||2012| ). This 
motivates a continued effort to attempt to confirm or deny 
a highly bursty star formation history for Sculptor. 



5.2.2 The star formation history of a population of dwarf 
galaxies 

An alternative approach to studying individual galaxies is to 
study a population statistically. If we assume that all dwarf 
galaxies of a given mass and type are statistically equiva- 
lent, then we can treat these has individual sample points 
along the SFH. If the SFH is smooth and continuous then 
the variance in the measured SFH for the population will 
be small; if, however, it is bursty then the variance will be 



larger. [Weisz et aL ( 2012 ) have recently used this idea to test 
a simple bursty star formation history against a population 
of 185 galaxies from the Spitzer Local Volume Legacy survey. 
They find that their more massive galaxies (M* ^ 1O 7 M0) 
are better fit (on average) by a smooth and continuous 
SFH, while the lower mass systems (M* ^ 1O 7 M0) favour 
a bursty SFH with bursts of ~ 10's of Myrs, inter-burst pe- 
riods of ~ 250 Myrs, and peak to trough burst amplitude 
ratios of ~ 30. These numbers are in total agreement with 
the SFH we obtained for our model with a cusp-core trans- 
formation (see Figure [3]). WLM is again interesting here. It 
has a stellar mass of 1.1 x 10 7 Mq a nd a dynamical mass 



of - 1O 1O M ( |Leaman et aL||2012) . As we have shown, 



it is comparable to our simulated halo, and statistically in 
the stellar mass range where bursty star formation in dwarf 
galaxies is observationally favored. 

It seems that all of the currently available data favor a 
bursty SFH for dwarf galaxies with stellar mass < 1O 7 M 
and dynamical mass < 1O 1O M . As we have shown here, 
such a bursty SFH will naturally heat the stars producing a 
hot stellar distribution even in isolated systems, consistent 
with recent observations of the isolated dwarf galaxy WLM. 
It also drives cusp-core transformations in the dark matter 
that can explain the now long-standing cusp-core problem 
( Flores fe Primack||1994| |Moore|1994| . 



6 DISCUSSION 

In this paper, we have simulated an isolated dwarf galaxy 
using the now traditional set-up of the cooling halo. Using 
a new implementation of stellar feedback within the RAM- 
SES code, we have modeled the formation of a dwarf galaxy 
of 10 10 M© halo mass, whose evolution has been shown to 
be strongly dominated by feedback processes. We have ob- 
served the formation of a 800 pc dark matter core, the size 
of which corresponds roughly to 40 resolution elements. We 
have analyzed in details the stellar distribution and kinemat- 
ics of our model galaxy, observable properties that compare 
successfully to the local isolated dwarf WLM. In our fidu- 
cial model, the total stellar mass is too large by one order of 
magnitude, suggesting that we need to lower the long-lived 



star formation efficiency, while keeping the same overall feed- 
back efficiency to match all constraints. We have shown that 
this might be indeed achieved, if one uses a very low value 
for the Schmidt law parameter e* c± 0.2%, together with a 
top-heavy initial mass function. 

Our simulation with feedback and AMR clearly con- 
firms previous works performed with various SPH codes, 
both for isolated haloes and in a cosmological context, at 
a somewhat lower resolution than the one used here. The 
temporal variations and systematic evolution of the inner 
slope of the dark matter profile correlate nicely with bursts 
seen in the SF history of our dwarf galaxy. We noted that 
there are several observational diagnostics which point to 
the cusp-flattening process driven by bursty star formation. 
The most obvious place to look is in star formation histories. 
Unfortunately for individual cases, it is unlikely to be possi- 
ble to achieve sufficient time resolution (e.g. Dolphin 2012; 
|de Boer et al.|[2012| . On the other hand, statistical studies 
of populations encouragingly point to bursts of exactly the 
right nature for cusp-flattening, and in agreement with our 
work, especially at low masses (Weisz et al. 2012). There 



may be less direct evidence of a cusp-flattening process in 
the final kinematics of the stars. Since stars behave as colli- 
sionless particles, like dark matter, they should be heated by 



the irreversible process proposed by Read & Gilmore 
|Mashchenko et ah] ( |2006j ) and |Pontzen fc Governato 



(2005), 



(2012) 



Indeed when cusp-flattening occurs, we also obtain hotter, 
more oblate stellar distributions with v/a ~ 1. More di- 
rectly related to the proposed mechanism, we also obtain a 
core in the stellar distribution. Promisingly, this also agrees 
well with the recent observations of the isolated dwarf galaxy 
WLM ( |Leaman et al.|2012 ). 

Although these observations tend to favor bursty SF 
histories as a plausible origin to cusp-core transformations, 
we would like to stress that this is only a necessary condi- 
tion for these to occur, not a sufficient one: other physical 
mechanisms could be at play to explain cusp-core transfor- 
mations, not related to baryonic processes, but, for example, 
to warm or self-interacting dark matter particles. 
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